/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | foam-extend: Open Source CFD
   \\    /   O peration     | Version:     4.1
    \\  /    A nd           | Web:         http://www.foam-extend.org
     \\/     M anipulation  | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
	This file is part of foam-extend.

	foam-extend is free software: you can redistribute it and/or modify it
	under the terms of the GNU General Public License as published by the
	Free Software Foundation, either version 3 of the License, or (at your
	option) any later version.

	foam-extend is distributed in the hope that it will be useful, but
	WITHOUT ANY WARRANTY; without even the implied warranty of
	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
	General Public License for more details.

	You should have received a copy of the GNU General Public License
	along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.

Description

\*---------------------------------------------------------------------------*/

#include "orientedSurface.H"
#include "triSurfaceTools.H"
#include "treeBoundBox.H"

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

defineTypeNameAndDebug(Foam::orientedSurface, 0);


// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //

// Return true if face uses edge from start to end.
bool Foam::orientedSurface::edgeOrder
(
	const labelledTri& f,
	const edge& e
)
{
	if
	(
		(f[0] == e[0] && f[1] == e[1])
	 || (f[1] == e[0] && f[2] == e[1])
	 || (f[2] == e[0] && f[0] == e[1])
	)
	{
		return true;
	}
	else
	{
		return false;
	}
}


// Return true if edge is used in opposite order in faces
bool Foam::orientedSurface::consistentEdge
(
	const edge& e,
	const labelledTri& f0,
	const labelledTri& f1
)
{
	return edgeOrder(f0, e) ^ edgeOrder(f1, e);
}


Foam::labelList Foam::orientedSurface::faceToEdge
(
	const triSurface& s,
	const labelList& changedFaces
)
{
	labelList changedEdges(3*changedFaces.size());
	label changedI = 0;

	forAll(changedFaces, i)
	{
		const labelList& fEdges = s.faceEdges()[changedFaces[i]];

		forAll(fEdges, j)
		{
			changedEdges[changedI++] = fEdges[j];
		}
	}
	changedEdges.setSize(changedI);

	return changedEdges;
}


Foam::labelList Foam::orientedSurface::edgeToFace
(
	const triSurface& s,
	const labelList& changedEdges,
	labelList& flip
)
{
	labelList changedFaces(2*changedEdges.size());
	label changedI = 0;

	// 1.6.x merge: using local faces.  Reconsider
	// Rewrite uses cached local faces for efficiency
	// HJ, 24/Aug/2010
	const List<labelledTri> lf =  s.localFaces();

	forAll(changedEdges, i)
	{
		label edgeI = changedEdges[i];

		const labelList& eFaces = s.edgeFaces()[edgeI];

		if (eFaces.size() < 2)
		{
			// Do nothing, faces was already visited.
		}
		else if (eFaces.size() == 2)
		{
			label face0 = eFaces[0];
			label face1 = eFaces[1];

			const labelledTri& f0 = lf[face0];
			const labelledTri& f1 = lf[face1];

			// Old.  HJ, 24/Aug/2010
//            const labelledTri& f0 = s[face0];
//            const labelledTri& f1 = s[face1];

			if (flip[face0] == UNVISITED)
			{
				if (flip[face1] == UNVISITED)
				{
					FatalErrorIn("orientedSurface::edgeToFace") << "Problem"
					    << abort(FatalError);
				}
				else
				{
					// Face1 has a flip state, face0 hasn't
					if (consistentEdge(s.edges()[edgeI], f0, f1))
					{
					    // Take over flip status
					    flip[face0] = (flip[face1] == FLIP ? FLIP : NOFLIP);
					}
					else
					{
					    // Invert
					    flip[face0] = (flip[face1] == FLIP ? NOFLIP : FLIP);
					}
					changedFaces[changedI++] = face0;
				}
			}
			else
			{
				if (flip[face1] == UNVISITED)
				{
					// Face0 has a flip state, face1 hasn't
					if (consistentEdge(s.edges()[edgeI], f0, f1))
					{
					    flip[face1] = (flip[face0] == FLIP ? FLIP : NOFLIP);
					}
					else
					{
					    flip[face1] = (flip[face0] == FLIP ? NOFLIP : FLIP);
					}
					changedFaces[changedI++] = face1;
				}
			}
		}
		else
		{
			// Multiply connected. Do what?
		}
	}
	changedFaces.setSize(changedI);

	return changedFaces;
}


void Foam::orientedSurface::walkSurface
(
	const triSurface& s,
	const label startFaceI,
	labelList& flipState
)
{
	// List of faces that were changed in the last iteration.
	labelList changedFaces(1, startFaceI);
	// List of edges that were changed in the last iteration.
	labelList changedEdges;

	while(true)
	{
		changedEdges = faceToEdge(s, changedFaces);

		if (debug)
		{
			Pout<< "From changedFaces:" << changedFaces.size()
				<< " to changedEdges:" << changedEdges.size()
				<< endl;
		}

		if (changedEdges.empty())
		{
			break;
		}

		changedFaces = edgeToFace(s, changedEdges, flipState);

		if (debug)
		{
			Pout<< "From changedEdges:" << changedEdges.size()
				<< " to changedFaces:" << changedFaces.size()
				<< endl;
		}

		if (changedFaces.empty())
		{
			break;
		}
	}
}


void Foam::orientedSurface::propagateOrientation
(
	const triSurface& s,
	const point& samplePoint,
	const bool orientOutside,
	const label nearestFaceI,
	const point& nearestPt,
	labelList& flipState
)
{
	//
	// Determine orientation to normal on nearest face
	//
	triSurfaceTools::sideType side = triSurfaceTools::surfaceSide
	(
		s,
		samplePoint,
		nearestFaceI,
		nearestPt,
		10*SMALL
	);

	if (side == triSurfaceTools::UNKNOWN)
	{
		// Non-closed surface. Do what? For now behave as if no flipping
		// nessecary
		flipState[nearestFaceI] = NOFLIP;
	}
	else if ((side == triSurfaceTools::OUTSIDE) == orientOutside)
	{
		// outside & orientOutside or inside & !orientOutside
		// Normals on surface pointing correctly. No need to flip normals
		flipState[nearestFaceI] = NOFLIP;
	}
	else
	{
		// Need to flip normals.
		flipState[nearestFaceI] = FLIP;
	}

	if (debug)
	{
		vector n = triSurfaceTools::surfaceNormal(s, nearestFaceI, nearestPt);

		Pout<< "orientedSurface::propagateOrientation : starting face"
			<< " orientation:" << nl
			<< "     for samplePoint:" << samplePoint << nl
			<< "     starting from point:" << nearestPt << nl
			<< "     on face:" << nearestFaceI << nl
			<< "     with normal:" << n << nl
			<< "     decided side:" << label(side)
			<< endl;
	}

	// Walk the surface from nearestFaceI, changing the flipstate.
	walkSurface(s, nearestFaceI, flipState);
}


bool Foam::orientedSurface::flipSurface
(
	triSurface& s,
	const labelList& flipState
)
{
	bool hasFlipped = false;

	// Flip tris in s
	forAll(flipState, faceI)
	{
		if (flipState[faceI] == UNVISITED)
		{
			FatalErrorIn
			(
				"orientSurface(const point&, const label, const point&)"
			)   << "unvisited face " << faceI
				<< abort(FatalError);
		}
		else if (flipState[faceI] == FLIP)
		{
			labelledTri& tri = s[faceI];

			label tmp = tri[0];

			tri[0] = tri[1];
			tri[1] = tmp;

			hasFlipped = true;
		}
	}
	// Recalculate normals
	s.clearOut();

	return hasFlipped;
}


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

// Null constructor
Foam::orientedSurface::orientedSurface()
:
	triSurface()
{}


// Construct from surface and point which defines outside
Foam::orientedSurface::orientedSurface
(
	const triSurface& surf,
	const point& samplePoint,
	const bool orientOutside
)
:
	triSurface(surf)
{
	orient(*this, samplePoint, orientOutside);
}


// Construct from surface. Calculate outside point.
Foam::orientedSurface::orientedSurface
(
	const triSurface& surf,
	const bool orientOutside
)
:
	triSurface(surf)
{
	// BoundBox calculation without localPoints
	treeBoundBox bb(surf.points(), surf.meshPoints());

	point outsidePoint = bb.max() + bb.span();

	orient(*this, outsidePoint, orientOutside);
}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

bool Foam::orientedSurface::orient
(
	triSurface& s,
	const point& samplePoint,
	const bool orientOutside
)
{
	bool anyFlipped = false;

	// Do initial flipping to make triangles consistent. Otherwise if the
	// nearest is e.g. on an edge inbetween inconsistent triangles it might
	// make the wrong decision.
	if (s.size() > 0)
	{
		// Whether face has to be flipped.
		//      UNVISITED: unvisited
		//      NOFLIP: no need to flip
		//      FLIP: need to flip
		labelList flipState(s.size(), UNVISITED);

		flipState[0] = NOFLIP;
		walkSurface(s, 0, flipState);

		anyFlipped = flipSurface(s, flipState);
	}


	// Whether face has to be flipped.
	//      UNVISITED: unvisited
	//      NOFLIP: no need to flip
	//      FLIP: need to flip
	labelList flipState(s.size(), UNVISITED);


	while (true)
	{
		// Linear search for nearest unvisited point on surface.

		scalar minDist = GREAT;
		point minPoint;
		label minFaceI = -1;

		forAll(s, faceI)
		{
			if (flipState[faceI] == UNVISITED)
			{
				const labelledTri& f = s[faceI];

				pointHit curHit =
					triPointRef
					(
					    s.points()[f[0]],
					    s.points()[f[1]],
					    s.points()[f[2]]
					).nearestPoint(samplePoint);

				if (curHit.distance() < minDist)
				{
					minDist = curHit.distance();
					minPoint = curHit.rawPoint();
					minFaceI = faceI;
				}
			}
		}

		// Did we find anything?
		if (minFaceI == -1)
		{
			break;
		}

		// From this nearest face see if needs to be flipped and then
		// go outwards.
		propagateOrientation
		(
			s,
			samplePoint,
			orientOutside,
			minFaceI,
			minPoint,
			flipState
		);
	}

	// Now finally flip triangles according to flipState.
	bool geomFlipped = flipSurface(s, flipState);

	return anyFlipped || geomFlipped;
}


// ************************************************************************* //
